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ABSTRACT 


A  set  of  basic  equations  necessary  to  model  the  acceleration  of  the  plasma 
produced  by  a  millipound  pulsed  plasma  thruster  is  presented  and  discussed.  Various 
simplifying  assumptions  and  approximations  typically  made  In  the  derivation  of 
mognetogasdynamic  equations  are  examined  in  light  of  available  experimental  data.  In 
particular,  it  is  shown  that  the  mean  free  path  and  Debye  shielding  length  are 
appropriate  to  the  assumption  of  a  neutral  continuum  with  pressure  and  Lorentz  force 
dominating  the  acceleration.  The  neglect  of  dispalcement  current  is  also  justified  and  a 
case  is  made  for  the  assumption  of  scalar  rather  than  tensor  thermal  and  electrical 
conductivities  and  the  neglect  of  viscous  effects. 

The  possibility  of  describing  the  flow  with  a  model  of  fewer  than  three  dimensions 
is  discussed  and  a  drastically  simplified  one-dimensiondl  model  is  presented,  along  with 
a  discussion  of  its  limitations  and  shortcomings. 


A  finite  difference  algorithm  for  the  one-dimensional  model  is  presented,  but 
difficulties,  apparently  arising  from  incorrect  boundary  conditions  prevented 
meaningful  numerical  experimentation. 
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I.  INTRODUCTION 


The  Teflon  PPT  has  been  around  for  some  fifteen  years.  Smoll  (approximately  one 
micropound  thrust)  versions  have  been  proved  feasible  on  satellites.  Larger  versions  are 
currently  being  used,  with  one  to  be  monitored  for  contamination  on  a  satellite  to  be 
launched  in  the  fall  of  1982.  The  thruster  of  primary  concern  in  this  report  is  the  much 
larger  (one  millipound  thrust)  model  which  is  still  in  the  process  of  being  developed. 

The  principles  of  PPT  operation  are  described  in  several  of  the  references  (see 
refs  I,  2)  and  only  a  very  brief  outline  is  given  here  for  completeness.  Bosicolly  a  PPT 
consists  of  a  pair  of  bar  electrodes  connected  by  a  capacitor  bank  with  external  power 
supply,  and  a  feed  system  which  forces  the  teflon  "fuel"  into  position  between  the 
electrodes.  This  is  shown  schematically  in  Figure  I.  The  capacitor  is  charged  between 
firings,  and  the  directly  coupled  electrodes  are  charged  os  well.  A  minute  amount  of 
material  is  introduced  into  the  electrode  gap  allowing  an  arc  discharge.  The  high 
temperature  arc  across  the  face  of  the  teflon  ablates,  dissociates  and  ionizes  a  smoll 
quantity  of  teflon  (producing  a  carbon-fluorine  plasma)  which  supports  the  continuing 
discharge.  The  current  flow  induces  a  high  magnetic  field  and  reacts  with  that  field  to 
accelerate  the  plasma  to  velocities  of  several  thousand  meters  per  second,  thus 
providing  the  thrust. 

The  very  high  exhaust  velocity  produces  one  of  the  PPT's  most  promising  features, 
a  specific  impulse  several  times  larger  than  that  of  chemical  rockets.  The  exhaust 
itself  poses  problems,  however.  A  serious  concern  of  the  satellite  designers  is  the 
contamination  potential  of  the  exhaust  plume.  The  plume  of  the  PPT  is  known  to 
contain  ionized  carbon  and  flourine  primarily,  but  additional  small  quantities  of  heavy 
molecules  (metals  and  metal  oxides  stripped  from  the  electrodes  ond  insulators  os  well 


as  flouroeorbans)  probably  also  exist.  It  is  important  to  determine  how  much  of  what 
exhaust  products  may  reach  sensitive  satellite  surfoces. 

The  obvious  potential  of  the  PPT  for  long  life,  dependable,  accurate,  and  very 
efficient  attitude  control  and  station  keeping  functions,  and  its  apparent  potential  as  a 
contamination  source,  have  been  the  subjects  of  a  considerable  amount  of  research. 
Reference  I  describes  much  of  the  work  pertinent  to  the  present  study,  and  that 
material  will  not  be  reviewed  in  detail  here,  except  to  note  that  references  2  and  3 
present  attempts  to  model  the  internal  and  external  plasma  flow  respectively.  These 
models  were  adopted  as  the  basis  of  the  investigation  presented  here,  and  will  be 
discussed  in  some  detail  below.  The  internal  model  was  developed  for  a  thruster  with 
an  entirely  different  feed  system  and  so  is  inapplicable  to  the  present  device  in  some 
respects. 

Important  new  experimental  data  are  now  available  from  investigations  performed 
at  the  Jet  Propulsion  Laboratory  and  the  Arnold  Engineering  Development  Center 
(references  4,  5  and  6)  which  have  enabled  additional  critical  evaluation  of  the  models. 
This  information  is  discussed  below  in  conjunction  with  its  consequences  with  regard  to 
various  aspect  of  the  model. 

2.  THRUSTER  OPERATION 


Figure  I  shows  the  millipound  thruster.  The  teflon  bars  ore  forced  in  against 
stops  machined  on  the  copper  anode.  The  igniter  is  fired,  releasing  a  minute  quantity  of 
ions  and  electrons  into  the  space  between  the  charged  electrodes  which  then  begi  i  to 
discharge.  The  discharge  initiates  the  ablation,  dissociation  and  ionization  of  the  face 
of  the  teflon  bars.  The  plasma  thus  formed  supports  the  continuing  current  discharge 
which  interacts  with  the  induced  magnetic  field  to  occeierate  the  plasma. 
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The  complete  electrical  circuit  (consisting  of  capacitor  bank,  electrodes,  and 
plasmoid,  plus  connections)  behaves  as  a  series  R,  L,  C  circuit.  That  is,  the  current 
exhibits  damped  sinusoidal  variation  as  shown  in  Figure  2.  The  contribution  of  the 
plasma  to  the  total  circuit  resistance  and  inductance  should  be  variable  since  the 
amount  of  plasmo  in  the  interelectrode  gap  changes  with  time.  (The  contribution  to  the 
capacitance  is  negligible.)  As  shown  in  the  figure,  however,  the  data  ogree  remarkably 
well  with  the  simple  damped  sine  wove  obtained  by  assuming  constant  circuit 
inductance  and  resistance.  The  coefficients  shown  in  the  figure  were  determined  by 
matching  the  data  to  the  equation  and  these  can  be  used  to  determine  the  circuit 
parameters  (see  reference  6) 

C  =  3.2(IO-8)farad 

R  =  1.3(10-2)  ohm 

L  =  9.7(10-8)  Henry 

It  is  possible  that,  at  least  during  the  first  cycle  (during  which  nearly  all  the  plasma  is 
formed),  the  plasma  components  of  resistance  and  inductance  are  small  compared  to 
those  of  the  fixed  portion  of  the  circuit,  so  that  the  totals  are  essentially  constont. 
During  the  second  current  cycle  very  little  additional  teflon  is  ablated,  ond  even  less  is 
ionized;  thus  the  resistance  is  significantly  increased  ond  the  damping  is  strongly 

t 

enhanced.  Similar  data  token  with  smaller  capacitors  installed  show  this  latter  effect 
more  strongly. 

If  one  assumes  that  the  amount  of  teflon  ablated  per  half-cycle  of  current 
discharge  is  proportional  to  the  integral  of  the  square  of  the  current  over  that  half- 
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cycle  (which  is  related  to  the  energy  consumed)  one  finds  that  some  92.1%  is  ablated 
during  the  first  half-cycle,  and  7.3%  in  the  next  half-cycle,  leaving  less  than  1%  for  the 
remainder  of  the  discharge.  As  will  be  discussed  in  Section  4.1,  below,  such  a  drastic 
difference  in  plasma  densities  makes  it  unlikely  that  a  single  theoretical  model  can  be 
developed  to  describe  the  entire  process. 

Further  support  for  this  assertion  was  sought  by  analyzing  the  (unpublished)  data 
obtained  by  discharging  the  thruster  through  an  aluminum  shorting  bar  across  the 
electrodes.  This  test  wos  run  in  the  atmosphere,  and  the  current  history  is  ogain  closely 
approximated  by  a  damped  sine  wave.  The  parameters  calculated  for  this  circuit  (with 
the  shorting  bar  replacing  the  plasma)  are 

C  =  3.2  (10-8)  farad 

R  =  IA  (10-2)  ohm 

L  =  2.1  (1 0~2)  Henry 

The  fact  that  the  resistance  agrees  well  with  thot  obtained  from  the  plasma 
forming  discharge  strengthens  the  assertion  that  the  resistance  of  the  plasma  itself  is 
small.  The  significant  (factor  of  two)  increase  in  the  inductonce  may  be  due  primarily 
to  the  increase  in  physical  size  of  the  current  loop.  (The  shorting  bar  was  located  at 
the  very  end  of  the  electrodes).  t 

The  field  vectors  in  Figure  I  indicate  the  direction  of  the  major  field  components 
during  the  first  holf-cycle  of  the  discharge.  It  should  be  noted  that  the  z-component  of 
the  magnetic  field  vector  changes  sign  somewhere  between  the  origin  ond  the  leading 


edge  of  the  plasmoid.  In  the  second  half -cycle  the  J  and  B  vectors  reverse,  but  the 
acceleration  vector  does  not.  Reference  6  provides  field  data  for  the  millipound 
thruster  which  appear  to  agree  qualitatively  with  those  for  the  micropound  thruster 
given  in  Reference  2. 


3.  MODEL  DEVELOPMENT  -  OVERVIEW 


A  complete  theoretical  model  of  the  PPT  operation  and  plume  can  be  partitioned 
into  an  internal  model  which  would  describe  the  circuit  discharge,  teflon  chemistry  and 
plasma  acceleration,  and  an  external  model  which  would  describe  the  plume  expansion 
beyond  the  ends  of  the  electrodes.  The  internal  model  must  realistically  represent  the 
physics  of  the  plasma  flow  using  appropriate  field  equations  for  the  fluid  flow 
parameters  (densities,  velocities,  pressures,  temperatures)  and  the  electromagnetic 
field  variables  (electric  field,  magnetic  field,  current  density,  etc.).  Also  required  are 
constitutive  equations  which  describe  the  nature  of  the  medium  (transport  coefficients, 
etc.),  equations  to  model  the  formation  of  the  plasmoid  (ablation,  ionization),  and 
boundary  conditions.  The  external  model  will  also  deal  with  plasma  flow.  However, 
since  the  plasmoid  expands  rapidly  into  a  vacuum,  the  density  quickly  becomes  so  low 
that  the  continuum  equations  of  the  internal  model  are  no  longer  appropriate.  Thus,  a 
different  approach  is  needed.  (This  may  also  be  true  for  the  interelectrode  flow  of  all, 
but  the  first  one  or  two  half-cycles  of  the  discharge.)  In  addition,  if  the  model  is  to  be 
used  in  conjunction  with  ground  tests,  chamber  effects  must  be  considered. 


4.  MODEL  DEVELOPMENT  -  INTERNAL  MODEL 


4.1  Field  Equations 
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Basically,  the  plasma  flow  variables  are  governed  by  the  time  dependent 
conservation  equations  (mass,  momentum,  energy)  and  the  Maxwell  electromagnetic 
field  equations.  Before  a  detailed  discussion  of  these,  however,  some  attention  should 
be  paid  to  the  medium  (i.e.,  the  nature  of  the  plasmoid). 

The  plasma  produced  by  the  PPT  discharge  is  a  complex  medium  which,  although 
the  subject  of  several  investigations,  has  yet  to  be  fully  characterized  or  understood. 
The  principal  constituents  are  carbon,  flourine,  and  flourocarbons  which  originate  in  the 
teflon.  These  atoms  and  molecules  may  be  neutral  or  ionized.  In  addition,  elements 
stripped  from  the  electrodes  and  other  parts  of  the  thruster  may  contribute. 

Reference  6  presents  spectroscopic  data  (in  the  UV  and  VUV  range)  and  identifies 
at  least  the  following  emitters  in  the  plasmoid  at  the  exit  plane. 

C®,  r  t  c»,  Fo,  F->,  CoS  Al+ 

Previous  studies  have  reported  additional  ions.  The  copper  is  apparently  from  the 
electrodes  while  the  aluminum  probably  originates  as  an  oxide  in  the  ceramic  insulator 
which  protects  the  sides  of  the  teflon  bars  from  the  hot  plasma. 

Several  observations  can  be  made  from  the  spectroscopic  data.  First,  the 
radiation  is  extremely  intense,  especially  in  the  VUV  range  where  it  is  estimated  to  be 
on  the  order  of  10^  times  the  solar  VUV  flux.  An  estimate  pf  the  temperature  is  also 
given  to  be  obout  40,000  K. 

The  spectra  exhibit  sharp,  well  defined  peaks  characte  istic  of  atomic  rather  than 
molecluar  emission  although  flourocarbons  do  emit  in  the  2300-3300  A  range  (reference 
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6).  This  would  seem  to  indicate  that  the  plasma  environment  ts  too  severe  for  the 
heavy  molecules  to  survive,  or  that  dissociation  is  complete. 

The  study  also  presents  data  obtained  by  monitoring  the  temporal  history  of 
several  individual  emission  lines  at  the  thruster  exit  plane.  All  of  the  lines  monitored 
exhibit  a  strong  pulse  between  eleven  and  nineteen  microseconds  into  the  cycle  and  a 
second  much  smaller  (by  as  much  as  a  factor  of  ten),  pulse  between  thirty  and  thiry  five 
microseconds.  The  species  so  monitored  include  aluminum,  neutral  and  ionized  carbon, 
and  flourine.  All  pulses  are  quite  similar  in  character  and  timing. 

Several  deductions  about  the  nature  of  the  plasma  can  be  made  from  the  above 
observations.  The  most  important  is  that  the  flow  is  dominated  by  collision  effects. 
That  this  must  be  the  case  is  evidenced  by  the  fact  that  atoms  with  vastly  different 
(even  zero)  charge  to  mass  ratio  pass  the  exit  plane  at  the  same  time,  and  apparently 
with  quite  similar  velocities.  In  addition,  this  information  implies  that,  althouqh  the 
plasmoid  is  highly  complex  mixture  of  species,  a  model  which  considers  it  to  be  a  three, 
or  even  one,  component  fluid  may  be  adequate  to  give  reasonable  accuracy. 

A  three  component  model  might  consist  of  neutrals,  ions,  and  electrons,  as 
suggested  in  reference  2.  There  it  was  proposed  that  on  average  molecular  weight  of 
16.67  be  used  for  the  ions  and  neutrals  (this  being  the  weighted  average  for  two 
flourines  to  one  carbon),  and  that  the  number  densities  of  ions  and  elections  are  equal 
(assuming  charge  neutrality  and  single  ionization). 

In  this  author's  opinion,  we  know  both  too  much  and  too  little  to  adopt  this  model. 
First,  it  is  apparent  that  single  ionization  is  not  a  qood  assumption.  Also,  it  is  quite 
likely  that  the  average  molecular  weights  of  the  ion  gas  and  the  neutral  gas  ore  not 
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equal,  since  the  ionization  potentials  of  carbon  and  f  lourlne  ore  quite  different.  On  the 
other  hand,  sufficient  data  are  not  available  to  enable  one  to  predict  with  any 
confidence  the  energy  transfer  between  the  various  gas  components.  Reference  2 
presents  much  pertinent  data,  such  as  dissociation  and  ionization  energies,  as  well  os  a 
model  for  ionization  rates,  etc.  Of  importance  here  is  the  fact  that  the  total  energy 
required  to  ablate,  dissociate  and  ionize  the  total  mass  of  teflon  is  a  small  fraction  of 
that  available,  or  that  required  to  accelerate  the  plasma  to  such  high  velocities.  Since 
reaction  rates,  collision  cross  sections,  etc.,  are  not  known,  it  seems  inappropriate  to 
include  such  refinements  in  a  model  at  this  time.  In  foct,  no  data  are  available  with 
which  to  compare  the  output  of  such  a  model.  It  thus  would  seem  preferable,  as  well  as 
much  simpler,  to  model  the  plasma  as  a  single  conducting  gas.  In  this  case,  the 
conservation  equations  consist  of  one  mass,  one  to  three  momentum  (depending  on 
dimensionality  of  the  model),  and  one  energy  equotion. 

There  are,  of  course,  inaccuracies  built  in  to  such  a  one  fluid  model.  For 
example,  doppler  shift  spectroscopy  measurements  (references  6  and  7)  seem  to 
indicate  somewhat  different  velocities  for  different  species.  These  data,  however,  are 
not  entirely  consistent  and  are  difficult  to  interpret.  Another  anomaly  is  demonstrated 
by  the  data  taken  at  different  positions  relative  to  the  thruster  axis  (Reference  5). 
These  show  considerable  variation  in  the  intensity  of  certain  emission  lines  with 
position.  Again,  interpretation  is  difficult. 

Some  approximate  calculations  have  been  performed  to  check  the  validity  of  the 
assertions  above.  The  mean  free  path  in  the  plasma  is  essentially  proportional  to  the 
inverse  of  the  number  density,  n,  whereas  the  Debye  length  is  proportional  to  the 
square  root  of  temperature/number  density.  Total  collision  cross  sections  of  carbon  and 
or  flourine  atoms  are  not  available  but  assumed  diameters  on  the  order  of  l(H0  m  give 


a  mean  free  path  X  of  about  10^/n  m.  The  Debye  length  Is  given  by 


V  «o*<T/nje2  1/2 


where  «0  is  the  free  space  permittivity,  K  is  Boltzman’s  constant,  T  is  temperature,  n, 
is  the  ion  number  density  and  e  is  the  electron  charge.  This  formula  yields. 

Xft=  69(T/n)l/2  meters 
with  T  in  Kelvin. 


To  approximate  the  number  density  in  the  ptasmoid  one  can  refer  to  photographic 
evidence.  At  II  sec  after  initiation  the  plasmoid  just  about  fills  the  interelectrode 
volume,  and  about  half  of  the  first  half-cycle  energy  has  been  deposited.  Assuming, 
then,  that  about  .7mg  of  material  with  an  overage  molecular  weight  of  16.67  occupies  a 
volume  of  1.5(10"**)  m3  yields  approximately  20023)  for  the  number  density  of  the  first 
plasmoid.  These  approximations  then  yield  a  mean  free  path  X  of  500-5)  m.  The 
additional  assumptions  that  the  plasmoid  is  about  20%  ionized  yields  o  Debye  length 
of  8(10-8)  m.  The  first  figure  alows  the  use  of  continuum  equations  as  X  is  small 
compared  to  physical  dimensions,  and  since  X  the  assumption  of  charge  neutraility 
is  also  valid. 

* 

It  must  be  noted,  however,  that  the  above  computations  are  volid  only  for  the 
first  plasmoid  (that  produced  by  the  first  half-cycle  of  the  discharge).  As  the  ablated 
mass  appears  to  decrease  by  an  order  of  magnitude  with  each  successive  half-cycle  the 
mean  free  path  will  increase  accordingly,  and  certainly  by  the  third  holf-cycle  the 
continuum  hypothesis  breaks  down.  The  Debye  length  also  increases,  but  inversely  as 
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V 


the  root  of  the  density,  Thus  a  model  developed  to  describe  the  first  plosmoid  may  not 
be  occurate  for  the  second  one  and  certainly  will  not  apply  to  any  subsequent  emissions. 
Unfortunately,  from  the  viewpoint  of  contamination  onalysis,  it  moy  well  be  that  the 
latter  portion  of  the  thruster  exhaust  is  of  prominent  importance. 

For  the  reasons  given  above  It  is  suggested  that  the  one  fluid  model  be  developed 
further.  For  such  a  model  the  governing  equations  respresent  conservation  of  mass, 
energy,  and  three  components  of  momentum,  coupled  with  Maxwell's  equations  of 
electromagnetics.  These  can  be  written  as: 

Pt  ♦  div  (p  v)  =  0 

e*  ♦  v  •  grad  e  =  conduction  ♦  dissipation  effects 

pvt  ♦  PY  *qrad  v  ♦  grad  p  =  Lorentz  bodv  force  ♦  viscous  forces 

div  E  sO 

Curl  E  =  -Bf 

div  B  =  O 

curl  B  =  p  (J  -  «0  Et>  . 

Here  p  -  gas  density 

v  =  gas  mean  velocity 

t  »  time  and  (  )t  represents  partial  differentiation  w.r.t.t 

n 


e  «  energy  per  unit  moss 


p  »  gas  pressure 

E  *  electric  field  intensity 

B  s  magnetic  flux  density 

J  s  current  density 

M  3  magnetic  permeability 

«o  3  vacuum  permittivity. 

To  these  equations  must  be  coupled  constitutive  equations  which  further  deserve  the 
medium  itself.  These  are  discussed  in  the  next  section. 

4.2  Constitutive  Equations 

In  the  most  general  one  fluid  plasma  models  with  stronq  magnetic  fields  present 
the  constitutive  equations  are  tensor  relationships.  That  is,  such  plasmas  are 
non- isotropic,  with  the  direction  of  the  mognetic  flux  density  vector  B  being  a 
"preferred  direction'’  in  the  sense  that  materiol  properties  irt  this  direction  differ  from 
those  in  perpendicular  directions.  The  properties  involved  include  both  electric  and 
thermal  conductivities  and  the  viscosity.  These  control  the  right  hand  sides  of  the 
momentum  and  energy  equations  above  and,  therefore,  some  preliminary  analysis  is  in 
order  at  this  time. 


Any  computation  of  the  plasma  transport  coefficients  is,  at  best,  an  estimate, 
but,  since  no  experimental  data  exist  for  a  plasma  resembling  that  produced  by  the 
PPT,  one  is  required  to  make  such  an  estimate.  Reference  8  presents  a  theory  for 
computing  transport  coefficients  for  a  two  component  plasmo  (singly  charged  ions  and 
electrons).  In  particular,  this  theory  allows  one  to  estimate  the  different  components 
of  the  transport  tensors. 

The  strength  of  the  magnetic  field  is  conveniently  described  by  the  cyclotron 
frequency  for  the  ions 

Oj  =  eB/mi 

where  e  is  the  ionic  charge  and  mj  the  mass.  This  is  the  frequency  at  which  a  single 
free  ion  would  orbit.  For  a  singly  charged  ion  of  atomic  mass  number  16.67  (mass  = 

2.88(10*26)  kg)  and  a  field  of  about  .5  Telsa  (which  is  appropriate  for  the  first  plasmoid) 
one  obtains 

0 1  *  3  (10^)  sec-1 

A  measure  of  the  collision  effect  is  a  collision  time  which  con  be  estimated  from 

,  ,  <0...m.v‘(KT)*/* 

*  rUC*  log  A 

0 

where  A  *  l.2H(lOr)  T*Un"** 

For  a  temperature  of  40,000  °K  and  a  number  density  of  1023  the  obove  yields  o 
collision  time  of  about  10*?  sec.  Since  the  orbital  frequency  (lt  *  3(10$)  sec"*,  it  is 


apparent  that  an  ion  would  experience  many  collisions  per  orbit.  This  calculation 
substantiates  the  assertion  above  that  the  process  is  collision  dominated,  and  also 
implies  that  a  scalar  approximation  for  the  transport  coefficients  is  appropriate.  Since 
the  collisions  are  randominizing  events,  the  "preference"  for  the  direction  parallel  to 
the  magnetic  field  is  lost.  When  an  ion  begins  to  react  to  the  presence  of  the  field  it 
suffers  o  collision  and  is  redirected. 

It  still  remains  to  provide  some  model  for  the  magnitudes  of  the  scalar  transport 
coefficients.  Theoretical  estimates  based  on  kinetic  theory  are  available,  but  all  such 
theories  presume  plasmas  of  a  much  less  complex  composition  thon  the  one  produced  by 
the  PPT.  Thus  any  model  chosen  will  impart  a  serious  uncertainty  on  the  calculations, 
unless  it  is  found  that  the  coefficeints  are  of  such  a  magnitude  that  variations  do  not 
alter  the  flow  in  a  significant  fashion.  Otherwise,  some  experimental  verification  of 
the  models  will  be  needed. 

Reference  8  provides  some  formulae  for  transport  coefficients,  based  on  a  two- 
component  plasmas  (i.e.,  highly  ionized),  as  follows 

0  s  ,7H  D-fe-Tt  =  electrical  conductivity 

me 

*  *  ^  *  Nierm°l  conductivity 

f  *  ,5  nKTTj  s  viscosity  * 

where  re  is  the  collision  time  for  the  electrons  and  is  equal  to  (Me/M;)*^  r\.  Using 
these  formulae  with  numerical  values  appropriate  to  the  first  plosmoid  yields  the 
following  order  of  magnitude  estimates 


0  s  I0*»  mhos/m 


«  *  10  J/sec  °K 
n  =  10-8  N-S/m2 

To  cxomine  the  validity  of  these  estimates  one  must  rely  on  some  very  rough 
comparisons  with  available  data.  For  example,  in  the  case  of  conductivity,  if  it 
assumed  that  the  piasmoid  provides  a  20cm 2  cross  section  for  current  flow  between  the 
electrodes  which  are  about  8cm  apart,  the  above  conductivity  yields  a  totol  resistance 
of  about  4(I0~3)Q.  This  is  about  one-fourth  of  the  total  circuit  resistance  as 
determined  by  the  current  discharge  measurements  (see  section  2.0  above),  which  is 
probably  high,  but  not  completely  unreasonable  according  to  the  eorlier  assertion  that 
the  plosma  resistance  is  small. 

The  viscosity  above,  along  with  appropriate  values  of  velocity  (10,000  m/s)  and 
density  (3XIO-3kg/m2)  gives  a  Reynolds  number 

Re  =  108 

which  implies  that  viscous  forces  ore  negligibly  small  except  within  an  extremely  thin 
boundary  layer  along  the  electrodes.  It  is  thus  proposed  to  ignore  viscous  terms  in  the 
model.  The  only  experimental  verification  of  this  assumptiog  is  that  photographs  of  the 
piasmoid  seem  to  indicate  foirly  uniform  velocity  profiles  which  would  not  occur  under 
the  action  of  strong  viscous  effects. 


With  nothing  better  then  o  rough  global  estimate  of  temperature,  no  similar 
conclusions  may  be  drown  concerning  heot  transfer.  In  fact,  considering  the  intensity 
of  both  the  UV  emission  and  visible  emission  from  the  plasmoid;  it  is  possible  that 
radiation  is  more  important  in  heot  transfer  than  thermal  conduction. 

Thus,  at  this  point,  the  choice  of  transport  coefficient  models  is  not  entirely 
clear.  A  reasonable  approach  would  seem  to  be  to  use  the  above  expressions  for 
electrical  and  thermal  conductivity  and  to  ignore  viscous  effects.  (Note  that,  unlike 
oridinary  gases,  plasmas  conduct  heat  primarily  through  the  electrons  and  therefore 
thermal  conduction  ond  viscous  diffusion  are  not  coupled  by  o  Prandtl  number  near 
unity). 

One  further  simplification  of  the  model  equations  is  usually  made  in 
mognetodydrodynamic  theory,  namely  the  neglect  of  the  displacement  current,  <0  Ef» 
compared  with  the  conduction  current  J.  The  magnitudes  of  the  variables  in  the  first 
PPT  plasmoid  are  completely  compatible  with  this  simplification  (which  is  quite 
significant),  but  some  caution  is  required.  In  the  vacuum  ahead  of  and  behind  the 
plasmoid  there  can  be  no  conduction  current.  The  electric  and  magnetic  fields, 
however,  do  not  vanish  (or  become  uniform).  Disturbances  propagate  with  the  velocity 
of  light,  which  is  much  higher  than  either  the  ion  velocity  or  the  Afven  wove 
propagation  in  the  plasma.  Thus,  the  MHD  equations  can  be  used  within  the  plasmoid, 
but  the  computation  domain  must  contain  only  thot  region  in  space  so  occupied.  The 
boundaries  of  the  computations  domain  thus  move  and  must  be  relocated  ot  each  time 
step.  The  boundary  conditions  applied  at  these  boundaries  must  consist  of  matching  the 
field  variables  across  the  plasmoid/vocuum  interface. 


38  Applying  the  simplifications  above  to  the  general  equations,  eliminating  the 

eleetrlc  field  using  the  generalized  Ohm’s  law  with  scoiar  conductivity 


E  =  B  x  V  ♦  J/e 

and  eliminating  the  current  density  using 
J  s  Xcurl  B 

(which  is  a  result  of  neglecting  displacement  current)  produces  a  familiar  set  of 
equations 

♦  div(pV)  s  0 

p vf  ♦  - grad  v  ♦  grad  v  grad  p  =  (E  ♦  V  ♦  B)  x  B 

pie  ♦  v  •  v/2)j  ♦  p  v  •  grad  (e  ♦  v  •  v/2)  v)  ♦  div  (p  v)  =  # (E  ♦  V  x  BVE  ♦  Txx 
To  these  must  be  added  the  algebraic  equations  of  state 

p  =  pRT 
e  *  CpT 

$ 

and  either  the  above  expressions  for  r  and  «  or  so.ne  other  model.  The  stote  equations 
here  assume  a  perfect  gas,  and  assuminq  a  monatomic  gas  with  atomic  weiqht  16.67 
gives  R  *  500,  Cp  s  2.5R. 
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4.3  One  Dimensional  Model 


This  model  should  provide  a  reasonable  approximation  to  reality.  For  the  purpose 
of  the  present  study,  the  model  was  simplified  still  further  for  numerical 
implementation.  First,  the  dimensionality  of  the  model  was  reduced  to  one  by 
neglecting  variations  in  the  y  and  z  directions  and  considering  only  the  x -component  of 
V,  the  y  -  component  of  J,  and  the  2 -component  of  B.  It  is  apparent  that,  while  these 
are  indeed  the  dominant  components  of  these  vectors,  the  one  dimensional  model  can 
provide  only  a  general  description  of  the  process.  Because  of  this  limitation  and  the 
uncertainty  in  the  transport  coefficients  and  ,  the  dissipation  (Joule  heating)  and 
heat  conduction  terms  were  dropped  and  temperature  was  taken  to  be  constant.  The 
one-dimensional,  isothermal  model  is  governed  by  the  equations 

Pt  *  (pu)*  =  S 
Mt*Px  =  j“BBx 
Bf  ♦  (uB)x  =  0 

p  *  pRT(T  constant) 

Here  the  source  term  S  must  be  Introduced  since  the  addition  of  plosma  due  to  the 
ablation  and  ionization  of  teflon  occurs  over  the  first*  few  centimeters  in  the 
x-direction.  In  the  full  three-dimensional  model  this  material  addition  would  appear  as 
a  boundary  condition  rather  than  in  the  model  equations. 
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The  source  function  is  one  of  the  weakest  points  in  the  model  development 
(whether  as  here  or  as  a  boundary  condition).  The  rate  of  ablation  of  the  teflon  must 
depend  on  the  energy  flux  ot  the  surface.  Unfortunately,  energy  transfer  in  the  plosmo 
is  not  well  understood.  Conductivity  models  have  not  been  verified  experimentally, 

temperatures  have  not  been  measured,  radiation  has  not  been  measured  in  the  IR  range, 
etc.  In  addition,  ablation  rates  are  known  only  at  much  lower  temperatures.  Thus,  the 
ablation  model  is  at  this  time  very  uncertain.  For  the  purpose  of  testing  the  one 
dimensional  isothermal  model,  it  was  assumed  that  the  source  term  in  the  mass 
conservation  equation  is  proportional  to  the  square  of  the  current  density  ot  a 
particular  value  of  x,  as  long  as  x<4.5cm  (the  extent  of  the  teflon  surface).  Even  here 
the  proportionality  factor  is  a  matter  of  conjecture  and  numerical  experimentation. 

It  should  be  noted  here  that  the  idea  of  coupling  the  oblation  rate  to  the  total 
current  flow  at  any  time  (as  proposed  in  reference  2)  is  not  appropriate,  especially  for 
the  side  feed  thruster  currently  under  consideration;  the  ablation  is  a  much  more  local 
phenomenon  as  the  teflon  fuel  cannot  "see"  the  bulk  of  the  plasma.  Thus,  even  if 
radiation  is  the  primary  heat  transfer  mechonism,  only  that  plasma  located  in  the 
narrow  "combustion  chamber"  between  the  teflon  bars  can  contribute,  while  after  a  few 
microseconds  the  bulk  of  the  current  flow  is  well  out  between  the  electrodes.  In  fact, 
according  to  the  current  flow  data  of  reference  6  the  current  flow  extends  well  past  the 
electrodes  as  the  plasmoid  emerges  from  the  thruster. 

The  oblation  model  suggested  here  effectively  ossumes«that  ablation,  dissociation, 
heating,  and  ionization  times  are  short  enouqh  to  be  neglected,  and  that  the  actual 
power  source  which  accomplishes  this  is  that  produced  by  the  current  flow.  In  a  two 
dimensional  flow  model  a  similar  device  could  be  used  (assuminq  uniformity  in  the  z 
direction)  but  in  a  full  three  dimensional  model  ablation  would  take  place  on  a  boundary 
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surface*  In  the  latter  case  additional  Investigation  into  exactly  what  portion  of  the 
plasma  contributes  to  the  local  heat  transfer  would  be  required. 

A  drastically  simplified  set  of  equations  has  been  proposed  to  describe  the  plasma 
motion  and  the  magnetic  field.  One  consequence  of  the  simplifications  made  can  be 
examined  without  recourse  to  numerical  implementation  of  the  model  (which  is 
impossible  without  a  set  of  initial  and  boundary  conditions).  The  neglect  of  the 
displacement  current  reduces  Faraday's  equation  to 

A  curl  B  =  J 

as  discussed  above,  and  allows  the  elimination  of  J  from  the  equations.  The  reduction 
of  the  system  to  the  one -dimensional  model  produces 

J*  -  J-Bx 

where  J  and  B  are  the  y  and  z  components  respectively.  In  fact,  although  symmetry 
through  the  z=0  plane  implies  that  B|  (the  x -component  of  B)  is  zero  it  cannot  be  true 
that  B|z  s  0  unless  the  medium  is  unbounded  in  the  z-direction  (as  the  one-dimensional 
model  implies).  Comparison  of  current  flow  and  mage  tic  field  data  from  reference  6 
indicates  that  the  current  density  is  not  proportional  to  the  magnetic  field  gradient 
along  the  axis  as  the  model  implies  (see  Figure  3).  This  could  be  result  of  the  failure  of 
the  one-dimensional  simplification  or  the  fact  that  the  probes  used  to  measure  the 
current  and  field  disturbed  the  flow,  it  is  not  believed  that  the  displacement  current 
could  be  responsible. 


4.4  Numerical  Model 


An  explicit,  forward  time  step,  finite  difference  olgorith  was  chosen  to 
numerically  test  the  simple  model  developed  above.  Although  a  characteristics  method 
or  on  implicit  finite  difference  method  miqht  be  more  appropriate  for  the  simple  model, 
these  would  be  very  difficult  to  generalize  to  the  multidimensional  dissipotive  model. 
The  specific  algorithm  chosen  was  a  Lax-Wendroff  type,  similor  to  that  presented  in 
reference  2.  The  equations  are  written  in  conservation  form 

wt  ♦  fx  =  0 

where  w  is  a  component  of  the  variables  vector  (  ,  u,  Ftt  and  the  f  components  ore 
functions  of  the  variables.  A  single  time  step  is  taken  in  two  parts,  and  central 
differences  are  used  for  the  x  derivatives.  Thus,  if  wjj  represents  the  varioble  w  at  the 
nth  grid  space  ar»d  jth  time  step,  the  equations  can  be  written: 

•Cl  •  ,/2(wJn*l  <fJn+l-fJn> 

Cl  *  <!/1> 

|J+1  yj  _  At  ,fJ*l/2  fj+l/2x 
n  =  n  Ax  *Tn+l/2  '  Vl/2' 

This  procedure  must,  of  course,  be  modified  for  the  density  equation  to  include  the 
source  term  in  which  B*  is  replaced  by  a  central  difference. 

Schemes  such  os  this  are  known  to  be  numerically  unstable  if  the  time  step  is 
longer  than  the  time  required  for  information  to  progate  over  one  grid  space.  Thus,  the 
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loco!  sound  speeds  ond  Alfven  wove  speeds  must  be  computed  ot  eoch  time  step,  ond  the 
next  At  chosen  accordingly.  Here  At  was  taken  to  be  be 


This  algorithm  produces  a  solution  in  only  o  triangular  sort  of  domain  in  the 
x,t-plane,  as  shown  in  Figure  k,  since  it  provides  no  value  for  wH  at  either  the  first  or 
last  x  grid  point.  These  values  must  be  determined  by  appropriate  boundary  conditions. 
The  analysis  of  boundary  conditions  was  not  a  port  of  this  project,  but  is  included  in  a 
subsequent  study  conducted  at  the  Air  Force  Rocket  Propulsion  Laboratory. 

A  computer  program  based  on  this  algorithm  and  mathematical  model  was  tested 
and  behaved  reasonably.  The  model  appears  to  be  stable,  but  low  densities  cause  very 
short  time  steps  (since  the  Alfven  speed  varies  inversely  with  density).  The  results 
produced  depend  strongly  on  the  initial  spatial  variations  assumed,  as  expected,  and  all 
that  can  be  concluded  here  is  that  initial  conditions  which  seem  to  be  compatible  with 
the  experimental  data  for  a  particular  time  give  rise  to  later  conditions  which  are 
altered  in  the  correct  direction. 

Some  attempt  to  impose  boundary  conditions  was  made,  but  without  success.  The 
boundary  conditions  chosen  caused  very  steep  gradients  of  the  variables  ot  the 
boundaries  and  these  did  trigger  numerico!  instabilities,  even  if  the  time  steps  were 

0 

shortened.  The  proper  handling  of  the  boundary  conditions  is  obviously  crucial  to  the 
model's  application,  and  these  conditions,  especially  those  at  the  moving  plcsma  front, 
ore  non-triviol. 


5.  CONCLUSIONS 


The  multidimensionol  model  presented  above,  with  scalar  conductivities  and  zero 
viscosity,  is  believed  to  provide  a  reasonable  description  of  the  acceleration  of  the  first 
plasmoid  ejected  by  a  firing  of  the  PPT.  The  same  model  may  not,  however,  be 
accurate  for  the  second  or  later  plasmoids.  The  one -dimensional  isothermal  model 
which  was  implemented  numerically  should  provide  a  crude  approximation  if  boundary 
and  initial  conditions  are  handled  properly,  and  this  should  be  used  in  additional  studies 
of  these  boundary  conditions. 

Considerable  additional  study  will  be  required  before  a  truly  predictive  internal 
model,  capable  of  providing  input  for  a  plume  model,  is  obtained.  The  work  presented 
here  is  merely  a  small  step  in  the  process.  The  present  model,  even  in  the  one¬ 
dimensional  form,  may,  however,  prove  valuable  for  thruster  performance  studies,  sirce 
nearly  all  of  the  thrust  must  be  produced  by  the  first  plasmoid. 
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Figure  4.  Solution  Domain  of  Numerical  Algorithm 


